home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / zroots.i < prev   
Text File  |  1995-07-26  |  4KB  |  127 lines
  1. /*
  2.    ZROOTS.I
  3.    Laguerre's method for finding roots of complex polynomial.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func zroots(a, imsort=)
  11. /* DOCUMENT zroots(a)
  12.      returns the vector of (complex) roots of the (complex) 
  13.      polynomial: Sum(a(i)*x^(i-1)) using Laguerre's method.
  14.      The roots are sorted in order of increasing real parts.  Call with
  15.      non-zero imsort keyword to sort into increasing imaginary parts.
  16.      See Numerical Recipes (Press, et. al., Cambridge Univ. Press 1988),
  17.      section 9.5.
  18.  */
  19. {
  20.   m= numberof(a)-1;   /* degree of poly */
  21.   EPS= 1.e-6;
  22.   EPS*= 2.*EPS;
  23.   is_complex= (structof(a)==complex);
  24.   if (!is_complex) a= a+0i;
  25.   ad= a;              /* copy of coeffs for successive deflation */
  26.   roots= array(0i, m);
  27.  
  28.   for (j=m ; j>=1 ; j--) {
  29.     /* Loop over each root to be found */
  30.     x= laguerre(ad,0i);     /* start at zero to favor smallest rem root */
  31.     if (abs(x.im)<=EPS*abs(x.re)) x= x.re+0i;
  32.     roots(j)= x;
  33.     b= ad(j+1);             /* forward deflation */
  34.     for (jj=j ; jj>=1 ; jj--) {
  35.       c= ad(jj);
  36.       ad(jj)= b;
  37.       b= x*b+c;
  38.     }
  39.     ad= ad(1:j);            /* poly has one lower degree */
  40.   }
  41.  
  42.   /* polish */
  43.   for (j=1 ; j<=m ; j++) roots(j)= laguerre(a,roots(j));
  44.  
  45.   /* sort roots */
  46.   re= roots.re;
  47.   im= roots.im;
  48.   if (!is_complex &&
  49.       allof(abs(im)<=EPS*abs(re))) return re(sort(re));
  50.  
  51.   /* sorting on multiple keys is difficult because all fast sorting
  52.      algorithms (such as the quicksort used by sort) randomize the
  53.      order of equal elements in the list to be sorted -- the following
  54.      is a modificatio n of the general deterministic sorting algorithm
  55.      implemented in msort.i, modified to ensure that conjugate roots
  56.      are recognized even when their real parts are not quite equal */
  57.   if (imsort) {
  58.     rank= re;
  59.     re= im;
  60.     im= rank;
  61.   }
  62.   list= rank= sort(re);
  63.   re= re(list);
  64.   rank(list)= (re(dif)>EPS*abs(re(zcen)))(cum);  /* see msort.i */
  65.   list= rerank= sort(im);
  66.   rerank(list)= indgen(m);
  67.   return roots(sort(rank+double(rerank)/m));
  68. }
  69.  
  70. func laguerre(a,x)
  71. /* DOCUMENT laguerre(a,x)
  72.      Given the coefficients a(1:m+1) of the m'th degree
  73.      complex polynomial Sum(a(i)*x^(i-1)) and a guess x, returns a root.
  74.      See Numerical Recipes (Press, et. al., Cambridge Univ. Press 1988),
  75.      section 9.5.
  76.  */
  77. {
  78.   EPSS=3.e-15;
  79.   MR=8;
  80.   MT=10;
  81.   MAXIT=MT*MR;
  82.   frac=[.5,.25,.75,.13,.38,.62,.88,1.]; /* Fractions to break limit cycles */
  83.  
  84.   m= numberof(a)-1;       /* degree of poly */
  85.   a= a+0i;                /* make complex */
  86.   if (m<2) return -a(1)/a(2);
  87.  
  88.   for (iter=1 ; iter<=MAXIT ; iter++) {
  89.     b= a(0);                    /* coefft of max degree */
  90.     d= f= 0i;
  91.     err= abs(b);
  92.     abx=abs(x);
  93.     for (j=m ; j>=1 ; j--) {    /* compute poly and 2 derivs */
  94.       f= f*x+d;
  95.       d= d*x+b;
  96.       b= b*x+a(j);
  97.       err= abs(b)+abx*err;
  98.     }
  99.     err*= EPSS;                 /* estimated round-off error in poly */
  100.     if (abs(b)<=err) return x;  /* we're on the root already */
  101.  
  102.     g= d/b;
  103.     g2= g*g;
  104.     h= g2-2.*f/b;
  105.     sq= sqrt((m-1)*(m*h-g2));
  106.     gp= g+sq;
  107.     gm= g-sq;
  108.     abp= abs(gp);
  109.     abm= abs(gm);
  110.     if (abp<abm) gp= gm;
  111.  
  112.     if (gp) dx= m/gp;
  113.     else dx= exp(log(1+abx)+1.i*iter);
  114.  
  115.     x1= x-dx;
  116.     if (abs(x-x1)<EPSS) return x1; /* converged to machine accuracy */
  117.  
  118.     if (iter%MT) x= x1;            /* next iteration */
  119.     else x= x-dx*frac(iter/MT);    /* break limit cycle? */
  120.   }
  121.   write, "WARNING: too many iterations in laguerre";
  122.   return x;
  123. }
  124.       
  125.  
  126.  
  127.